# ---------------------------------------------------------------- #
# The Journal of Politics
# When Do Citizens Resist The Use of AI Algorithms in Public Policy?
#
# Replication Code for Decision Type Experiment
# 
# Author: Shir Raviv
# Date: Feb 19, 2025
# ---------------------------------------------------------------- #
# Clear environment and free memory ---------------
rm(list = ls())
gc()

# Load required packages -----------
required_packages <- c("dplyr", "readr", "tidyr", "expss","scales",'ggeffects','lme4','forcats',
                       "vtable","xtable",'stargazer',"lme4",'texreg','broom','rlang',
                       "rio", "broom","ggplot2", 'tableone')
for (pkg in required_packages) {
  if (!require(pkg, character.only = TRUE)) {
    install.packages(pkg, dependencies = TRUE)
    library(pkg, character.only = TRUE)
  }
}

# Set working directory and define export paths ---------------
#setwd("insert here.. /JOP Replication")
plotpath <- "figures/"

# Load Survey Data  ---------------
df <- readRDS("DT_wide_data.RDS")
long_df <- readRDS("DT_long_data.RDS")

long_df_bwcomp<-long_df%>%
  filter(component_matrix=="Between")

long_df_withincomp<-long_df%>%
  filter(component_matrix=="Within")

# Table (A-1) Summary statistics for the survey sample ---------
# Create age summary
age_summary <- df %>%
  mutate(age_group = case_when(
    age_cat == "18 - 24" ~ "18-24",
    age_cat == "25 - 34" ~ "25-34",
    age_cat == "35 - 44" ~ "35-44",
    age_cat == "45 - 54" ~ "45-54",
    age_cat %in% c("55 - 64", "65 +") ~ "55+"
  )) %>%
  group_by(age_group) %>%
  summarise(
    Sample_N = n(),
    Sample_Percent = round(n()/nrow(df)*100, 1)
  ) %>%
  mutate(
    Category = "Age",
    Subcategory = age_group,
    Population_Percent = case_when(
      age_group == "18-24" ~ 13.0,
      age_group == "25-34" ~ 19.0,
      age_group == "35-44" ~ 22.0,
      age_group == "45-54" ~ 18.0,
      age_group == "55+" ~ 28.0
    )
  )

# Create race summary
race_summary <- df %>%
  mutate(race_group = case_when(
    race_cat == "White" ~ "White",
    race_cat == "Hispanic/Latino" ~ "Hispanic",
    race_cat == "Black/African" ~ "Black",
    race_cat == "Asian" ~ "Asian",
    TRUE ~ "Other"
  )) %>%
  filter(race_group != "Other") %>%
  group_by(race_group) %>%
  summarise(
    Sample_N = n(),
    Sample_Percent = round(n()/nrow(df)*100, 1)
  ) %>%
  mutate(
    Category = "Race",
    Subcategory = race_group,
    Population_Percent = case_when(
      race_group == "White" ~ 60.7,
      race_group == "Hispanic" ~ 18.1,
      race_group == "Black" ~ 13.4,
      race_group == "Asian" ~ 7.8
    )
  )

# Create gender summary
gender_summary <- df %>%
  group_by(female_fac) %>%
  summarise(
    Sample_N = n(),
    Sample_Percent = round(n()/nrow(df)*100, 1)
  ) %>%
  mutate(
    Category = "Gender",
    Subcategory = female_fac,
    Population_Percent = case_when(
      female_fac == "Male" ~ 48.0,
      female_fac == "Female" ~ 52.0
    )
  )

education_summary <- df %>%
  group_by(low_edu_somecol_cat) %>%
  summarise(
    Sample_N = n(),
    Sample_Percent = round(n()/nrow(df)*100, 1)
  ) %>%
  mutate(
    Category = "Education",
    Subcategory = low_edu_somecol_cat,
    Population_Percent = case_when(
      low_edu_somecol_cat == "Associate’s Degree or Higher" ~ 38.3,
      low_edu_somecol_cat == "Some College or Less" ~ 61.7)
  )

# Combine all summaries
summary_table <- bind_rows(
  age_summary,
  race_summary,
  gender_summary,
  education_summary
) %>%
  select(Category, Subcategory, Population_Percent, Sample_N, Sample_Percent)
summary_table

# Table (A-3) Summary of statistics of Perceived Accuracy, Fairness, and Appropriateness ---------
# Function to summarize data
summarize_data <- function(data, decision_col, outcome_col) {
  decision_sym <- sym(decision_col)
  outcome_sym <- sym(outcome_col)
  
  data %>%
    dplyr::select(!!decision_sym, !!outcome_sym) %>%
    filter(!is.na(!!outcome_sym), !is.na(!!decision_sym)) %>%
    group_by(!!decision_sym) %>%
    summarize(
      mean = mean(!!outcome_sym),
      n = n(),
      sd = sd(!!outcome_sym),
      se = sd / sqrt(n),
      .groups = 'drop'
    ) %>%
    mutate(Type = as.character(!!decision_sym))
}

# Applying function to different contexts and outcomes
data_app_mean_domain <- bind_rows(
  summarize_data(df, 'edu_decision_type', 'DTexp_app_mech_edu'),
  summarize_data(df, 'poli_decision_type', 'DTexp_app_mech_poli'),
  summarize_data(df, 'child_decision_type', 'bin_app_mech_child'),
  .id = "context"
)

data_fair_mean_domain <- bind_rows(
  summarize_data(df, 'edu_decision_type', 'DTexp_fair_mech_edu'),
  summarize_data(df, 'poli_decision_type', 'DTexp_fair_mech_poli'),
  summarize_data(df, 'child_decision_type', 'DTexp_fair_mech_child'),
  .id = "context"
)

data_accu_mean_domain <- bind_rows(
  summarize_data(df, 'edu_decision_type', 'DTexp_accur_mech_edu'),
  summarize_data(df, 'poli_decision_type', 'DTexp_accur_mech_poli'),
  summarize_data(df, 'child_decision_type', 'DTexp_accur_mech_child'),
  .id = "context"
)

# Combine all data into one dataframe
data_mean_domain <- bind_rows(
  list(Appropriateness = data_app_mean_domain,  Fairness = data_fair_mean_domain, Accuracy = data_accu_mean_domain),
  .id = "consideration"
)

# Adjuste labels for context (domain) and consideration
data_mean_domain$context <- factor(data_mean_domain$context, levels = c("1", "2", "3"), labels = c("Education", "Policing", "Child Welfare"))

data_mean_domain$consideration <- factor(data_mean_domain$consideration, levels = c("Appropriateness", "Fairness", "Accuracy"), labels = c("Appropriateness","Fairness","Accuracy" ))

data_mean_domain<-data_mean_domain%>%
  dplyr::select(consideration, context, Type, mean,n, sd, se)

# LaTeX table
latex_table <- xtable(data_mean_domain)
print(latex_table, type = "latex", include.rownames = FALSE)
 
# Table (A-4) Summary of statistics of binary outcomes ------
pairwise_bw_component <- long_df_bwcomp %>%
  group_by(T_decisiontype, T_domain) %>%
  summarize(
    mean_appropriate = mean(outcome_approp_bin3, na.rm = TRUE),
    se_appropriate = sd(outcome_approp_bin3, na.rm = TRUE) / sqrt(n()),
    mean_fairness = mean(outcome_fair_bin3, na.rm = TRUE),
    se_fairness = sd(outcome_fair_bin3, na.rm = TRUE) / sqrt(n()),
    mean_accuracy = mean(outcome_accur_bin3, na.rm = TRUE),
    se_accuracy = sd(outcome_accur_bin3, na.rm = TRUE) / sqrt(n())
  )
pairwise_bw_component

# Table (A-5) Balance tests --------

# Define the list of demographic variables
demographics <- c("female_fac", "age_cat", "low_edu_somecol_cat", "white", "high_diglit")

# Balance table for Policing
balance_policing <- df %>%
  filter(DTexp_order_policing == 1) %>%
  select(all_of(demographics), poli_decision_type)

# Create a balance table by treatment groups (poli_decision_type)
table_policing <- CreateTableOne(vars = demographics,
                                 strata = "poli_decision_type",
                                 data = balance_policing,
                                 factorVars = demographics)
print(table_policing, showAllLevels = TRUE)

# Balance table for Education
balance_education <- df %>%
  filter(DTexp_order_education == 1) %>%
  select(all_of(demographics), edu_decision_type)

table_education <- CreateTableOne(vars = demographics,
                                  strata = "edu_decision_type",
                                  data = balance_education,
                                  factorVars = demographics)
print(table_education, showAllLevels = TRUE)

# Balance table for Child Welfare
balance_childwelfare <- df %>%
  filter(DTexp_order_child == 1) %>%
  select(all_of(demographics), child_decision_type)

table_childwelfare <- CreateTableOne(vars = demographics,
                                     strata = "child_decision_type",
                                     data = balance_childwelfare,
                                     factorVars = demographics)

print(table_childwelfare, showAllLevels = TRUE)


# Table (A-6) Decision type and perceived appropriateness, by policy domain --------

# --- Domain: Education ---
# Base model for Education (subset where DTexp_order_education == 1)
lm_app_edu1_order1 <- lm(bin_app_mech_edu ~ edu_assisting + edu_individuals,
                         data = subset(df, DTexp_order_education == 1))

# Model with demographic controls and digital literacy
lm_app_edu2_order1 <- lm(bin_app_mech_edu ~ edu_assisting + edu_individuals +
                           age + female_fac + low_edu_somecol_cat + white + high_diglit,
                         data = subset(df, DTexp_order_education == 1))

# Model with demographic controls and prior AI knowledge
lm_app_edu3_order1 <- lm(bin_app_mech_edu ~ edu_assisting + edu_individuals +
                           age + female_fac + low_edu_somecol_cat + white + AI_have_priorknow,
                         data = subset(df, DTexp_order_education == 1))

# Rename coefficients for education models
names(lm_app_edu1_order1$coefficients) <- c('Constant', 'Assisting', 'Individuals')
names(lm_app_edu2_order1$coefficients) <- c('Constant', 'Assisting', 'Individuals',
                                            "Age", "Female", 
                                            "Some college or less", "White", "High tech literacy")
names(lm_app_edu3_order1$coefficients) <- c('Constant', 'Assisting', 'Individuals',
                                            "Age", "Female",
                                            "Some college or less", "White", "Prior Knowledge")


# --- Domain: Policing ---
# Base model for Policing (subset where DTexp_order_policing == 1)
lm_app_poli1_order1 <- lm(bin_app_mech_poli ~ poli_assisting + poli_individuals,
                          data = subset(df, DTexp_order_policing == 1))

# Model with demographic controls and digital literacy
lm_app_poli2_order1 <- lm(bin_app_mech_poli ~ poli_assisting + poli_individuals +
                            age + female_fac + low_edu_somecol_cat + white + high_diglit,
                          data = subset(df, DTexp_order_policing == 1))

# Model with demographic controls and prior AI knowledge
lm_app_poli3_order1 <- lm(bin_app_mech_poli ~ poli_assisting + poli_individuals +
                            age + female_fac + low_edu_somecol_cat + white + AI_have_priorknow,
                          data = subset(df, DTexp_order_policing == 1))

# Rename coefficients for policing models
names(lm_app_poli1_order1$coefficients) <- c('Constant', 'Assisting', 'Individuals')
names(lm_app_poli2_order1$coefficients) <- c('Constant', 'Assisting', 'Individuals',
                                             "Age", "Female",
                                             "Some college or less", "White", "High tech literacy")
names(lm_app_poli3_order1$coefficients) <- c('Constant', 'Assisting', 'Individuals', 
                                             "Age", "Female",
                                             "Some college or less", "White", "Prior Knowledge")


# --- Domain: Child Welfare ---
# Base model for Child Welfare (subset where DTexp_order_child == 1)
lm_app_child1_order1 <- lm(bin_app_mech_child ~ child_assisting + child_individuals,
                           data = subset(df, DTexp_order_child == 1))

# Model with demographic controls and digital literacy
lm_app_child2_order1 <- lm(bin_app_mech_child ~ child_assisting + child_individuals +
                             age + female_fac + low_edu_somecol_cat + white + high_diglit,
                           data = subset(df, DTexp_order_child == 1))

# Model with demographic controls and prior AI knowledge
lm_app_child3_order1 <- lm(bin_app_mech_child ~ child_assisting + child_individuals +
                             age + female_fac + low_edu_somecol_cat + white + AI_have_priorknow,
                           data = subset(df, DTexp_order_child == 1))

# Rename coefficients for child welfare models  
names(lm_app_child1_order1$coefficients) <- c('Constant', 'Assisting', 'Individuals')
names(lm_app_child2_order1$coefficients) <- c('Constant', 'Assisting', 'Individuals',
                                              "Age", "Female",
                                              "Some college or less", "White", "High tech literacy")
names(lm_app_child3_order1$coefficients) <- c('Constant', 'Assisting', 'Individuals',
                                              "Age", "Female", 
                                              "Some college or less", "White", "Prior Knowledge")


# --- 4. Create LaTeX Table using texreg ---
# Reorder the models to match the desired column headers:
# Education: Columns 1-3; Policing: Columns 4-6; Child Welfare: Columns 7-9.
models_list <- list(
  lm_app_edu1_order1,  # Column (1): Education, Base model
  lm_app_edu2_order1,  # Column (2): Education, Demog & Digital literacy
  lm_app_edu3_order1,  # Column (3): Education, Demog & Prior Knowledge
  lm_app_poli1_order1,  # Column (4): Policing, Base model
  lm_app_poli2_order1,  # Column (5): Policing, Demog & Digital literacy
  lm_app_poli3_order1,  # Column (6): Policing, Demog & Prior Knowledge
  lm_app_child1_order1, # Column (7): Child Welfare, Base model
  lm_app_child2_order1, # Column (8): Child Welfare, Demog & Digital literacy
  lm_app_child3_order1  # Column (9): Child Welfare, Demog & Prior Knowledge
)

# Specify custom model names to reflect domain and model specification.
custom_names <- c(
  "Education: Base",
  "Education: Demog & Digital",
  "Education: Demog & Prior",
  "Policing: Base",
  "Policing: Demog & Digital",
  "Policing: Demog & Prior",
  "Child Welfare: Base",
  "Child Welfare: Demog & Digital",
  "Child Welfare: Demog & Prior"
)

# LaTeX table 
texreg(models_list,
       digits = 3,
       custom.model.names = custom_names,
       caption = "Decision type and perceived appropriateness, by policy domain",
       label = "tab:dc_app_perceived",
       stars = c(0.1, 0.05, 0.01, 0.001),
       include.adjrs = FALSE,
       include.aic = FALSE,
       include.bic = FALSE,
       include.loglik = FALSE,
       include.deviance = FALSE,
       custom.gof.names = c("R$^2$","Observations"),
       booktabs = TRUE
)


# Table (A-7) Decision type and appropriateness, full sample controlling for order --------
# Primary Outcome Measure (Binary: 3 Last Categories) 
lm_app_edu_full1 <- lm(bin_app_mech_edu ~ edu_assisting + edu_individuals + DTexp_order_education +
                         age + female_fac + low_edu_somecol_cat + white + high_diglit, data = df)
lm_app_poli_full1 <- lm(bin_app_mech_poli ~ poli_assisting + poli_individuals + DTexp_order_policing +
                          age + female_fac + low_edu_somecol_cat + white + high_diglit, data = df)
lm_app_child_full1 <- lm(bin_app_mech_child ~ child_assisting + child_individuals + DTexp_order_child +
                           age + female_fac + low_edu_somecol_cat + white + high_diglit, data = df)

coef_names2 <- c('Constant', 'Assisting', 'Individuals', 'Order',
                 "Age", "Female", "Some college or less", "White", "Tech literacy")
names(lm_app_edu_full1$coefficients) <- coef_names2
names(lm_app_poli_full1$coefficients) <- coef_names2
names(lm_app_child_full1$coefficients) <- coef_names2

# Alternative Outcome Measure (Binary: 2 Last Categories) ---
library(dplyr)
df <- df %>%
  mutate(bin_app_mech_edu2 = case_when(DTexp_app_mech_edu %in% c(7, 6) ~ 1,
                                       DTexp_app_mech_edu < 6 ~ 0),
         bin_app_mech_poli2 = case_when(DTexp_app_mech_poli %in% c(7, 6) ~ 1,
                                        DTexp_app_mech_poli < 6 ~ 0),
         bin_app_mech_child2 = case_when(DTexp_app_mech_child %in% c(7, 6) ~ 1,
                                         DTexp_app_mech_child < 6 ~ 0))

lm_app_edu_full2 <- lm(bin_app_mech_edu2 ~ edu_assisting + edu_individuals + DTexp_order_education +
                         age + female_fac + low_edu_somecol_cat + white + high_diglit, data = df)
lm_app_poli_full2 <- lm(bin_app_mech_poli2 ~ poli_assisting + poli_individuals + DTexp_order_policing +
                          age + female_fac + low_edu_somecol_cat + white + high_diglit, data = df)
lm_app_child_full2 <- lm(bin_app_mech_child2 ~ child_assisting + child_individuals + DTexp_order_child +
                           age + female_fac + low_edu_somecol_cat + white + high_diglit, data = df)

names(lm_app_edu_full2$coefficients) <- coef_names2
names(lm_app_poli_full2$coefficients) <- coef_names2
names(lm_app_child_full2$coefficients) <- coef_names2

# Alternative Outcome Measure (Seven-Point Scale) ---
lm_app_edu_full3 <- lm(DTexp_app_mech_edu ~ edu_assisting + edu_individuals + DTexp_order_education +
                         age + female_fac + low_edu_somecol_cat + white + high_diglit, data = df)
lm_app_poli_full3 <- lm(DTexp_app_mech_poli ~ poli_assisting + poli_individuals + DTexp_order_policing +
                          age + female_fac + low_edu_somecol_cat + white + high_diglit, data = df)
lm_app_child_full3 <- lm(DTexp_app_mech_child ~ child_assisting + child_individuals + DTexp_order_child +
                           age + female_fac + low_edu_somecol_cat + white + high_diglit, data = df)

names(lm_app_edu_full3$coefficients) <- coef_names2
names(lm_app_poli_full3$coefficients) <- coef_names2
names(lm_app_child_full3$coefficients) <- coef_names2

# LaTeX Table  ---
# Education: Columns 1-3; Policing: Columns 4-6; Child welfare: Columns 7-9.
models_list <- list(
  lm_app_edu_full1,  # Column (1): Education, Binary (3 last cat)
  lm_app_edu_full2,  # Column (2): Education, Binary (2 last cat)
  lm_app_edu_full3,  # Column (3): Education, 7-pt scale
  lm_app_poli_full1,  # Column (4): Policing, Binary (3 last cat)
  lm_app_poli_full2,  # Column (5): Policing, Binary (2 last cat)
  lm_app_poli_full3,  # Column (6): Policing, 7-pt scale
  lm_app_child_full1, # Column (7): Child welfare, Binary (3 last cat)
  lm_app_child_full2, # Column (8): Child welfare, Binary (2 last cat)
  lm_app_child_full3  # Column (9): Child welfare, 7-pt scale
)

custom_names <- c(
  "Education: Binary (3 last cat)",
  "Education: Binary (2 last cat)",
  "Education: 7-pt scale",
  "Policing: Binary (3 last cat)",
  "Policing: Binary (2 last cat)",
  "Policing: 7-pt scale",
  "Child welfare: Binary (3 last cat)",
  "Child welfare: Binary (2 last cat)",
  "Child welfare: 7-pt scale"
)

texreg(models_list,
       digits = 3,
       custom.model.names = custom_names,
       caption = "Decision type and appropriateness, full sample controlling for order",
       label = "tab:dc_exp_full_order",
       stars = c(0.1, 0.05, 0.01, 0.001),
       include.adjrs = FALSE,
       include.aic = FALSE,
       include.bic = FALSE,
       include.loglik = FALSE,
       include.deviance = FALSE,
       custom.gof.names = c("R$^2$","Observations"),
       booktabs = TRUE
)


# Table (A-8) Decision type and appropriateness, controlling for inattentiveness --------

# Recode Inattentiveness Variables ---

# (a) Response-time attentiveness: convert 'inattentive' variable to a factor.
df$inattentive_time <- factor(df$inattentive)

# (b) Attention check: create a variable based on an embedded item in the matrix.
df$inattentive_matrix <- 0
df$inattentive_matrix[df$app_mech_attention_both != 5] <- 1
df$inattentive_matrix <- factor(df$inattentive_matrix)

# Estimate Models ---

# OLS Models controlling for response-time inattentiveness
lm_app_edu_full_attentive1 <- lm(bin_app_mech_edu ~ edu_assisting + edu_individuals + DTexp_order_education +
                                   age + female_fac + low_edu_somecol_cat + white + high_diglit +
                                   inattentive_time, data = df)
lm_app_poli_full_attentive1 <- lm(bin_app_mech_poli ~ poli_assisting + poli_individuals + DTexp_order_policing +
                                    age + female_fac + low_edu_somecol_cat + white + high_diglit +
                                    inattentive_time, data = df)
lm_app_child_full_attentive1 <- lm(bin_app_mech_child ~ child_assisting + child_individuals + DTexp_order_child +
                                     age + female_fac + low_edu_somecol_cat + white + high_diglit +
                                     inattentive_time, data = df)

# Rename coefficients to more meaningful names.
coef_names2 <- c('Constant', 'Assisting', 'Individuals', 'Order', "Age", "Female",
                 "Some college or less", "White", "Tech literacy", "Inattentive (time)")
names(lm_app_edu_full_attentive1$coefficients) <- coef_names2
names(lm_app_poli_full_attentive1$coefficients) <- coef_names2
names(lm_app_child_full_attentive1$coefficients) <- coef_names2

# OLS Models controlling for inattentiveness (Matrix item)
lm_app_edu_full_attentive2 <- lm(bin_app_mech_edu ~ edu_assisting + edu_individuals + DTexp_order_education +
                                   age + female_fac + low_edu_somecol_cat + white + high_diglit +
                                   inattentive_matrix, data = df)
lm_app_poli_full_attentive2 <- lm(bin_app_mech_poli ~ poli_assisting + poli_individuals + DTexp_order_policing +
                                    age + female_fac + low_edu_somecol_cat + white + high_diglit +
                                    inattentive_matrix, data = df)
lm_app_child_full_attentive2 <- lm(bin_app_mech_child ~ child_assisting + child_individuals + DTexp_order_child +
                                     age + female_fac + low_edu_somecol_cat + white + high_diglit +
                                     inattentive_matrix, data = df)

coef_names3 <- c('Constant', 'Assisting', 'Individuals', 'Order', "Age", "Female",
                 "Some college or less", "White", "Tech literacy", "Inattentive (matrix)")
names(lm_app_edu_full_attentive2$coefficients) <- coef_names3
names(lm_app_poli_full_attentive2$coefficients) <- coef_names3
names(lm_app_child_full_attentive2$coefficients) <- coef_names3

# Logistic Regression Models (without Inattentiveness controls)
glm_app_edu_logfull1 <- glm(bin_app_mech_edu ~ edu_assisting + edu_individuals + DTexp_order_education +
                              age + female_fac + low_edu_somecol_cat + white + high_diglit, 
                            data = df, family = binomial)
glm_app_poli_logfull1 <- glm(bin_app_mech_poli ~ poli_assisting + poli_individuals + DTexp_order_policing +
                               age + female_fac + low_edu_somecol_cat + white + high_diglit, 
                             data = df, family = binomial)
glm_app_child_logfull1 <- glm(bin_app_mech_child ~ child_assisting + child_individuals + DTexp_order_child +
                                age + female_fac + low_edu_somecol_cat + white + high_diglit, 
                              data = df, family = binomial)

coef_names4 <- c('Constant', 'Assisting', 'Individuals', 'Order',
                 'Age', 'Female', 'Some college or less', 'White', 'Tech literacy')
names(glm_app_edu_logfull1$coefficients) <- coef_names4
names(glm_app_poli_logfull1$coefficients) <- coef_names4
names(glm_app_child_logfull1$coefficients) <- coef_names4

# LaTeX Table ---

# Combine all models in desired order:
models_list <- list(
  lm_app_edu_full_attentive1,  # Education: OLS (time)
  lm_app_edu_full_attentive2,  # Education: OLS (matrix)
  glm_app_edu_logfull1,        # Education: Logistic
  lm_app_poli_full_attentive1,  # Policing: OLS (time)
  lm_app_poli_full_attentive2,  # Policing: OLS (matrix)
  glm_app_poli_logfull1,        # Policing: Logistic
  lm_app_child_full_attentive1, # Child Welfare: OLS (time)
  lm_app_child_full_attentive2, # Child Welfare: OLS (matrix)
  glm_app_child_logfull1        # Child Welfare: Logistic
)

# Specify custom model names for clarity:
custom_names <- c(
  "Education: OLS (time)",
  "Education: OLS (matrix)",
  "Education: Logistic",
  "Policing: OLS (time)",
  "Policing: OLS (matrix)",
  "Policing: Logistic",
  "Child Welfare: OLS (time)",
  "Child Welfare: OLS (matrix)",
  "Child Welfare: Logistic"
)

# Generate the LaTeX table:
texreg(models_list,
       digits = 3,  # show three decimals
       custom.model.names = custom_names,
       caption = "Decision type and appropriateness, controlling for inattentiveness",
       label = "tab:dc_exp_bw_inattentive_logistic",
       stars = c(0.1, 0.05, 0.01, 0.001),
       include.adjrs = FALSE,
       include.aic = TRUE,
       include.bic = FALSE,
       include.loglik = TRUE,
       include.deviance = FALSE,
       # Only include Observations and R^2 in the table
       custom.gof.names = c("R$^2$", "Observations", "Log Likelihood", "Akaike Inf. Crit."),
       booktabs = TRUE
)

# Figure (3) Effects of decision type on perceived appropriateness of ADS, across domains ------------
# Model Estimation by Domain ---
# Education model (subset where DTexp_order_education == 1)
lm_app_edu1 <- lm(bin_app_mech_edu ~ edu_assisting + edu_individuals,                      
                  data = subset(df, DTexp_order_education == 1))
# Policing model (subset where DTexp_order_policing == 1)
lm_app_poli1 <- lm(bin_app_mech_poli ~ poli_assisting + poli_individuals,
                   data = subset(df, DTexp_order_policing == 1))
# Child Welfare model (subset where DTexp_order_child == 1)
lm_app_child1 <- lm(bin_app_mech_child ~ child_assisting + child_individuals,
                    data = subset(df, DTexp_order_child == 1))

names(lm_app_edu1$coefficients) <- c('Constant','Objective: \n Assisting \n (ref: Sanctioning)', 'Target: \n Individuals \n (ref: Collectives)')
names(lm_app_poli1$coefficients) <- c('Constant','Objective: \n Assisting \n (ref: Sanctioning)', 'Target: \n Individuals \n (ref: Collectives)')
names(lm_app_child1$coefficients) <- c('Constant','Objective: \n Assisting \n (ref: Sanctioning)', 'Target: \n Individuals \n (ref: Collectives)')

# Tidy model outputs and combine results ---
m.lm_app_edu <-  tidy(lm_app_edu1, conf.int = TRUE)
m.lm_app_poli <-  tidy(lm_app_poli1, conf.int = TRUE)
m.lm_app_child <-  tidy(lm_app_child1, conf.int = TRUE)

app_domains_type <- bind_rows(
  m.lm_app_poli,
  m.lm_app_child,
  m.lm_app_edu,
  .id="domain")
app_domains_type$domain<-factor(app_domains_type$domain, 
                                levels=c(1,2,3), 
                                labels=c("Policing","Child Welfare","Education"))

app_domains_type<-subset(app_domains_type, 
                         term!="Constant" ) 
app_domains_type$term<-factor(app_domains_type$term)
app_domains_type

# Plot ---
fig3<-
  ggplot(app_domains_type, 
         aes(term, estimate, fill=term, colour=term, 
             label=round(estimate, digits = 2)))+ 
  scale_fill_manual(values = c("#3A6958", "#84A188"))+
  scale_colour_manual(values = c("#3A6958", "#8B181C"))+
  geom_hline(yintercept=0, linetype=2, lwd=0.9, size=2, 
             colour = "#ACBBBF", alpha=0.9) +
  geom_errorbar(stat = "identity", alpha = 0.5, 
                position = position_dodge(width = 0.15),
                aes(ymin=estimate - 1.96*std.error, 
                    ymax=estimate + 1.96*std.error),
                lwd=1.8, width = 0)+
  geom_errorbar(stat = "identity", alpha = 0.9, 
                position = position_dodge(width = 0.15),
                aes(ymin=estimate - 1.64*std.error, 
                    ymax=estimate + 1.64*std.error), 
                lwd=1.5, width=0) +
  geom_point(stat = "identity", alpha = 0.9,  size = 3, 
             position = position_dodge(width = 0.15)) + coord_flip() + 
  facet_grid(~domain, scale="free_y")+
  theme_minimal()+theme(legend.position = "none", 
                        panel.background = element_rect(fill = "white", colour = "gray30"), 
                        text = element_text(family = "Times New Roman", size=15.5),
                        strip.text = element_text(family = "Times New Roman", size=16),
                        axis.text = element_text(family = "Times New Roman"),
                        plot.background = element_rect(fill = "white", colour = "white"))+
  labs(x = "", y = " ", title="",  
       colour="", fill="", shape="", group="" ) 
fig3

#  Save the Figure ---
ggsave(fig3,
       filename = paste0(plotpath, "fig3.tif"),
       device = "tiff",
       height = 4, width = 9, dpi = 600)


ggsave(fig3,
       filename = paste0(plotpath, "fig3.png"),
       device = "png",
       height = 5, width = 9, dpi = 600)

# Figure (4) Perceived Fairness versus Accuracy, by Decision Type and Issue Area ----------
# Using data (long version) collected from the second component of the matrix 
pairwise_comparisons <- long_df_withincomp %>%
  group_by(T_decisiontype, T_decision) %>%
  summarize(
    mean_appropriate = mean(outcome_approp_bin3, na.rm = TRUE),
    mean_fairness = mean(outcome_fair_bin3, na.rm = TRUE),
    mean_accuracy = mean(outcome_accur_bin3, na.rm = TRUE),
    t_test = list(t.test(outcome_fair_bin3, outcome_accur_bin3, paired = TRUE)),
    p_value = sapply(t_test, function(x) x$p.value)
  )

pairwise_comparisons$T_decisiontype <- factor(pairwise_comparisons$T_decisiontype, 
                                              levels = c("assist collectives", "sanction collectives", 
                                                         "assist individuals", "sanction individuals"))

# Reshape the data to long format
pairwise_comparisons_long <- pairwise_comparisons %>%
  pivot_longer(cols = c(mean_fairness, mean_accuracy), names_to = "Measure", values_to = "Value")


pairwise_comparisons_long$T_decision <- factor(pairwise_comparisons_long$T_decision, levels = c(
  "Fire stations",  "Homeless Shelters",
  "Construction",  "Immigration",
  "Food stamps",  "Study assistance", 
  "Restraining order" ,  "Criminal Sentencing" # Ensure these two are ordered correctly
))

pairwise_comparisons_long$T_decisiontype <- factor(pairwise_comparisons_long$T_decisiontype, 
                                                   levels = c("assist collectives", "sanction collectives", 
                                                              "assist individuals", "sanction individuals"))

fairness_color <- "#8B181C"  # 
accuracy_color <- "#77A68C"  # 

fig4<-ggplot(pairwise_comparisons_long, aes(x = Value, y = reorder(T_decision, -as.numeric(T_decisiontype)), group = T_decision)) + 
  geom_line(aes(color = Measure), size = 0.5, colour="gray50") +  # Line connecting points
  geom_point(aes(color = Measure), size = 3.5) +
  labs(title = "",x = "", y = "") +  # No axis labels, as they're in the plot itself
  scale_color_manual(values = c(mean_fairness = fairness_color, 
                                mean_accuracy = accuracy_color),
                     labels = c("Accuracy", "Fairness"), name = "") +  # Empty legend title
  scale_x_continuous(
    breaks = seq(0.2, 0.6, by = 0.1), labels = percent_format(accuracy = 1),
    limits = c(0.2, 0.6))+ theme_minimal() +
  theme_minimal()+theme(legend.position = "bottom", 
                        panel.background = element_rect(fill = "white", colour = "gray30"), 
                        text = element_text(family = "Times New Roman", size=14.5),
                        strip.text = element_text(family = "Times New Roman", size=12),
                        axis.text = element_text(family = "Times New Roman"),
                        plot.background = element_rect(fill = "white", colour = "white"))+
  facet_grid(T_decisiontype ~ ., scales = "free_y", space = "free_y", 
             labeller = labeller(T_decisiontype = c(
               "assist collectives" = "Assist \n Collectives",
               "sanction collectives" = "Sanction \n Collectives",
               "assist individuals" = "Assist \n Individuals",
               "sanction individuals" = "Sanction \n Individuals" ))) 

# save plots
fig4
ggsave(fig4,
       filename = paste0(plotpath, "fig4.png"),
       device = "png",
       height = 5, width = 7, dpi = 600)


ggsave(fig4,
       filename = paste0(plotpath, "fig4.tif"),
       device = "tiff",
       height = 5, width = 7, dpi = 600)


# Figure (A-3) Decision type and predicted views on ADS, Between Subjects Component  -----
# Estimate mixed-effects models with random intercepts for domain and respondent
# Model for perceived appropriateness (binary outcome, 3 last categories)
model_lpm_app <- lmer(outcome_approp_bin3 ~ T_decisiontype + (1 | T_domain) + (1 | ResponseId),
                      data = long_df_bwcomp)

# Model for perceived fairness (binary outcome)
model_lpm_fair <- lmer(outcome_fair_bin3 ~ T_decisiontype + (1 | T_domain) + (1 | ResponseId),
                       data = long_df_bwcomp)

# Model for perceived accuracy (binary outcome)
model_lpm_accur <- lmer(outcome_accur_bin3 ~ T_decisiontype + (1 | T_domain) + (1 | ResponseId),
                        data = long_df_bwcomp)


# Obtain predicted values for T_decisiontype from each model 
lm1_app_predict   <- as.data.frame(ggpredict(model_lpm_app, terms = "T_decisiontype"))
lm1_fair_predict  <- as.data.frame(ggpredict(model_lpm_fair, terms = "T_decisiontype"))
lm1_accur_predict <- as.data.frame(ggpredict(model_lpm_accur, terms = "T_decisiontype"))

# Combine predictions 
data_predict <- bind_rows(
  lm1_app_predict   %>% mutate(Consideration = "Perceived Appropriate"),
  lm1_fair_predict  %>% mutate(Consideration = "Perceived Fairness"),
  lm1_accur_predict %>% mutate(Consideration = "Perceived Accuracy")
)

data_predict$Consideration <- factor(data_predict$Consideration, 
                                     levels = c("Perceived Appropriate", 
                                                "Perceived Fairness", 
                                                "Perceived Accuracy"))

# plot

plotA3<-
  data_predict %>%
  ggplot()+aes(y=predicted,x=x,
               colour=Consideration,fill=Consideration, shape=Consideration)+
  scale_fill_manual(values = list("#8B181C", "#2F5D62", "#5E8B7E" ))+
  scale_color_manual(values = list("#8B181C", "#2F5D62", "#5E8B7E" ))+
  geom_point(size=3,  position = position_dodge(width = 0.15)) +
  geom_errorbar(aes(ymin = predicted-std.error*1.96, ymax = predicted+std.error*1.96), 
                alpha=0.4, width = 0, size=0.8,  position = position_dodge(width = 0.15)) +
  geom_errorbar(aes(ymin = predicted-std.error*1.64, ymax = predicted+std.error*1.64), 
                width = 0, size=0.8,  position = position_dodge(width = 0.15)) +
  labs(x = "Decision Type
       ", y = "Predicted Value
       ", colour="", fill="", shape="", group="" ) +
  theme_minimal() + theme(legend.position = "bottom")
plotA3

#Save figure
ggsave(plotA3,
       filename = paste0(plotpath, "figA3.tif"),
       device = "tiff",
       height = 4, width = 9, dpi = 600)


ggsave(plotA3,
       filename = paste0(plotpath, "figA3.png"),
       device = "png",
       height = 5, width = 9, dpi = 600)


# Table (A-9) Mixed-Effects Linear Regressions for Perceived Appropriateness -----------
# Convert variables to factors to ensure proper treatment in the models.
long_df_bwcomp$T_domain           <- factor(long_df_bwcomp$T_domain)
long_df_bwcomp$T_type_assisting   <- factor(long_df_bwcomp$T_type_assisting)
long_df_bwcomp$T_type_individuals <- factor(long_df_bwcomp$T_type_individuals)
long_df_bwcomp$BWorder            <- factor(long_df_bwcomp$BWorder)

# Model 1: Dependent variable = outcome_approp_bin3 (Binary: 3 Last Categories)
model_mix1 <- lmer(outcome_approp_bin3 ~ T_type_assisting + T_type_individuals + BWorder +
                     T_domain + birthyear + female_fac + low_edu_somecol_cat + white + high_diglit +
                     (1 | ResponseId), data = long_df_bwcomp)

# Model 2: Dependent variable = outcome_approp_bin2 (Binary: 2 Last Categories)
model_mix2 <- lmer(outcome_approp_bin2 ~ T_type_assisting + T_type_individuals + BWorder +
                     T_domain + birthyear + female_fac + low_edu_somecol_cat + white + high_diglit +
                     (1 | ResponseId), data = long_df_bwcomp)

# Model 3: Dependent variable = outcome_approp (Continuous 7-Point Scale)
model_mix3 <- lmer(outcome_approp ~ T_type_assisting + T_type_individuals + BWorder +
                     T_domain + birthyear + female_fac + low_edu_somecol_cat + white + high_diglit +
                     (1 | ResponseId), data = long_df_bwcomp)

# LaTeX Table
stargazer(model_mix1, model_mix2, model_mix3,
          type = "latex",
          title = "Mixed-Effects Linear Regressions for Perceived Appropriateness",
          label = "tab:mix_effects_approp",
          dep.var.labels = c("Binary Outcome: 3 Last Categories",
                             "Binary Outcome: 2 Last Categories",
                             "7-Point Scale Outcome"),
          dep.var.caption = "Dependent Variables",
          covariate.labels = c("Assisting (Treatment)",
                               "Individuals (Treatment)",
                               "Order: Condition 2",
                               "Order: Condition 3",
                               "Domain: Education",
                               "Domain: Policing",
                               "Birth Year",
                               "Female",
                               "Low Education (Some College or Less)",
                               "White",
                               "High Digital Literacy"),
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          star.char = c("†", "*", "**", "***"),
          notes = "†p<0.1; *p<0.05; **p<0.01; ***p<0.001",
          notes.append = FALSE,
          # Additional formatting options:
          column.labels = c("(1)", "(2)", "(3)"),
          model.numbers = FALSE,
          align = TRUE,
          no.space = TRUE)



# Table (A-10) Full results from Figure A-3 -------
# Convert key variables to factors
long_df_bwcomp$T_domain <- factor(long_df_bwcomp$T_domain)
long_df_bwcomp$BWorder   <- factor(long_df_bwcomp$BWorder)

# Estimate mixed-effects models with random intercepts for domain and respondent
# Model for perceived appropriateness (binary outcome, 3 last categories)
model_lpm_app <- lmer(outcome_approp_bin3 ~ T_decisiontype + (1 | T_domain) + (1 | ResponseId),
                      data = long_df_bwcomp)

# Model for perceived fairness (binary outcome)
model_lpm_fair <- lmer(outcome_fair_bin3 ~ T_decisiontype + (1 | T_domain) + (1 | ResponseId),
                       data = long_df_bwcomp)

# Model for perceived accuracy (binary outcome)
model_lpm_accur <- lmer(outcome_accur_bin3 ~ T_decisiontype + (1 | T_domain) + (1 | ResponseId),
                        data = long_df_bwcomp)

# Table ---
stargazer(model_lpm_app, model_lpm_fair, model_lpm_accur,
          type = "latex",
          title = "Mixed-Effects Models for Predicted Views on ADS",
          label = "tab:between_subjects_component",
          covariate.labels = c("Assist Individuals",
                               "Sanction Collectives",
                               "Sanction Individuals", "Constant"),
          omit.stat = c("adj.rsq", "ser", "f"),
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          star.char = c("†", "*", "**", "***"),
          notes = "†p<0.1; *p<0.05; **p<0.01; ***p<0.001",
          notes.append = FALSE)

# Table (A-11) Summary of statistics of Perceived Accuracy, Fairness, and Appropriateness ---------
#  Define a Function to Summarize Data  ---
# Compute descriptive summaries for perceived  Appropriateness ---
long_df_withincomp$T_decisiontype <- factor(long_df_withincomp$T_decisiontype,
                                            levels = c("assist individuals", "assist collectives", "sanction individuals", "sanction collectives"),
                                            labels = c("Assist Individuals", "Assist Collectives", "Sanctioning Individuals", "Sanctioning Collectives")
)
# Rename labels in the T_decision
long_df_withincomp$T_decision_new <- ifelse(long_df_withincomp$T_decision == "Construction", "Illegal Construction",
                                            ifelse(long_df_withincomp$T_decision == "Immigration", "Illegal Work",
                                                   long_df_withincomp$T_decision))  


approp_summary <- long_df_withincomp %>%
  group_by(T_decisiontype, T_decision_new) %>%
  summarise(
    mean = mean(outcome_approp, na.rm = TRUE),
    n = sum(!is.na(outcome_approp)),
    sd = sd(outcome_approp, na.rm = TRUE),
    se = sd / sqrt(n))

# Compute descriptive summaries for perceived  Fairness ---
fairness_summary <- long_df_withincomp %>%
  group_by(T_decisiontype, T_decision_new) %>%
  summarise(
    mean = mean(outcome_fair, na.rm = TRUE),
    n = sum(!is.na(outcome_fair)),
    sd = sd(outcome_fair, na.rm = TRUE),
    se = sd / sqrt(n))

# Compute descriptive summaries for perceived accuracy ---
accuracy_summary <- long_df_withincomp %>%
  group_by(T_decisiontype, T_decision_new) %>%
  summarise(
    mean = mean(outcome_accur, na.rm = TRUE),
    n = sum(!is.na(outcome_accur)),
    sd = sd(outcome_accur, na.rm = TRUE),
    se = sd / sqrt(n))

# Combine  ---

data_types_combined <- bind_rows(approp_summary, fairness_summary,accuracy_summary,
                                 .id = "consideration")
data_types_combined$consideration<-factor(data_types_combined$consideration, levels=c(1,2,3),
                                          labels=c("Appropriateness","Fairness","Accuracy"))

consideration_order <- c("Appropriateness", "Fairness", "Accuracy")

decision_order <- c(
  "Food stamps", "Study assistance",          # Assist Individuals
  "Fire stations", "Homeless Shelters",        # Assist Collectives
  "Criminal Sentencing", "Restraining order",  # Sanction Individuals
  "Illegal Construction", "Illegal Work"    # Sanction Collectives
)
data_types_combined$consideration <- factor(data_types_combined$consideration, levels = consideration_order, ordered = TRUE)
data_types_combined$T_decision_new <- factor(data_types_combined$T_decision_new, levels = decision_order, ordered = TRUE)

data_types_combined <- data_types_combined %>%
  arrange(consideration, T_decision_new)

# LaTex Table  ---
xtable(data_types_combined)



# Table (A-12) Summary of statistics of binary outcomes ---------
pairwise_comparisons <- long_df_withincomp %>%
  group_by(T_decisiontype, T_decision_new) %>%
  summarize(
    mean_appropriate = mean(outcome_approp_bin3, na.rm = TRUE),
    se_appropriate = sd(outcome_approp_bin3, na.rm = TRUE) / sqrt(n()),
    mean_fairness = mean(outcome_fair_bin3, na.rm = TRUE),
    se_fairness = sd(outcome_fair_bin3, na.rm = TRUE) / sqrt(n()),
    mean_accuracy = mean(outcome_accur_bin3, na.rm = TRUE),
    se_accuracy = sd(outcome_accur_bin3, na.rm = TRUE) / sqrt(n())
  )

pairwise_comparisons
# Table (A-13) Decision type and appropriateness, controlling for issue area (within-subjects component)-------------
long_df_withincomp$WTorder<-factor(long_df_withincomp$WTorder)
# OLS Models
model_lpm1 <- lm(outcome_approp_bin3 ~ T_type_individuals + T_type_assisting +T_type_individuals * T_type_assisting+
                   assist_ind_decision_type + 
                   assist_col_decision_type + 
                   sanc_ind_decision_type + 
                   sanc_col_decision_type + 
                   WTorder+ factor(ResponseId), data = long_df_withincomp)

model_lpm2 <- lm(outcome_approp_bin2 ~ T_type_individuals + T_type_assisting +T_type_individuals * T_type_assisting+
                   assist_ind_decision_type + 
                   assist_col_decision_type + 
                   sanc_ind_decision_type + 
                   sanc_col_decision_type + 
                   WTorder+ 
                   factor(ResponseId), data = long_df_withincomp)


model_lpm3 <- lm(outcome_approp ~ T_type_individuals + T_type_assisting +T_type_individuals * T_type_assisting+
                   assist_ind_decision_type + 
                   assist_col_decision_type + 
                   sanc_ind_decision_type + 
                   sanc_col_decision_type + 
                   WTorder+ factor(ResponseId), data = long_df_withincomp)

# Mixed-Effects Models
model1 <- lmer(outcome_approp_bin3 ~ T_type_individuals + T_type_assisting +T_type_individuals * T_type_assisting+
                 assist_ind_decision_type + 
                 assist_col_decision_type + 
                 sanc_ind_decision_type + 
                 sanc_col_decision_type + 
                 WTorder+ 
                 (1 | ResponseId), data = long_df_withincomp)

model2 <- lmer(outcome_approp_bin2 ~ T_type_individuals + T_type_assisting +T_type_individuals * T_type_assisting+
                 assist_ind_decision_type + 
                 assist_col_decision_type + 
                 sanc_ind_decision_type + 
                 sanc_col_decision_type + 
                 WTorder+ 
                 (1 | ResponseId), data = long_df_withincomp)


model3 <- lmer(outcome_approp ~ T_type_individuals + T_type_assisting +T_type_individuals * T_type_assisting+
                 assist_ind_decision_type + 
                 assist_col_decision_type + 
                 sanc_ind_decision_type + 
                 sanc_col_decision_type + 
                 WTorder+ 
                 (1 | ResponseId), data = long_df_withincomp)

# LaTeX Table ---

stargazer(model_lpm1,model_lpm2,model_lpm3,
          model1,model2,model3, 
          type="text", omit="factor\\(ResponseId\\)",
          dep.var.labels=c("Appropriate (3 cat)", "Appropriate (2 cat)", "Appropriateness (7 cat)",
                           "Appropriate (3 cat)", "Appropriate (2 cat)", "Appropriateness (7 cat)"),
          covariate.labels = c(
            "Individuals",
            "Assisting",
            "Issue area: Study assistance",
            "Issue area: Shelters for homeless",
            "Issue area: Sentence",
            "Issue area: Illegal work",
            "Order 2","Order 3","Order 4",
            "Individuals × Assisting"),
          omit.stat = c('adj.rsq', 'ser', 'f'),
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001), star.char = c('†', '*', '**', '***'),
          notes = '†p<0.1; *p<0.05; **p<0.01; ***p<0.001', notes.append = F)



# Table (A-14) Pairwise Comparisons of Fairness and Accuracy by decision type and issue area--------
# By type of decision on average 
pairwise_comparisons_within_type_average <- long_df_withincomp %>%
  group_by(T_decisiontype) %>%
  summarize(
    mean_appropriateness = mean(outcome_approp_bin3, na.rm = TRUE),
    mean_fairness = mean(outcome_fair_bin3, na.rm = TRUE),
    mean_accuracy = mean(outcome_accur_bin3, na.rm = TRUE),
    t_test = list(t.test(outcome_fair_bin3, outcome_accur_bin3, paired = TRUE)),
    p_value = sapply(t_test, function(x) x$p.value))

avg_table <- pairwise_comparisons_within_type_average %>%
  mutate(T_decision = "Average") %>%
  select(T_decisiontype, T_decision, mean_fairness, mean_accuracy, p_value)

# By type of decision issue
pairwise_comparisons_within_type_issue <- long_df_withincomp %>%
  group_by(T_decisiontype, T_decision) %>%
  summarize(
    mean_appropriateness = mean(outcome_approp_bin3, na.rm = TRUE),
    mean_fairness = mean(outcome_fair_bin3, na.rm = TRUE),
    mean_accuracy = mean(outcome_accur_bin3, na.rm = TRUE),
    t_test = list(t.test(outcome_fair_bin3, outcome_accur_bin3, paired = TRUE)),
    p_value = sapply(t_test, function(x) x$p.value))

issue_table <- pairwise_comparisons_within_type_issue %>%
  select(T_decisiontype, T_decision, mean_fairness, mean_accuracy, p_value)

# Combine the average and issue-level tables and format the numeric columns
combined_table <- bind_rows(avg_table, issue_table) %>%
  arrange(T_decisiontype, T_decision) %>%
  mutate(
    mean_fairness = sprintf("%.3f", mean_fairness),
    mean_accuracy = sprintf("%.3f", mean_accuracy),
    p_value = sprintf("%.3f", p_value)
  )

combined_table
